# !/bin/bash

# 主要是将indel和snp的vcf文件进行合并, 增加merge_vcf函数, 生成的结果文件id_samplename.merge.vcf.gz, 考虑到肿瘤数据库分析
# version = 4.0

usage(){
    echo "
Example:
ctdna.sh s=<sample>
"
}

default(){
    a3=GATCGGAAGAGCACACGTCT
    a5=AATGATACGGCGACCACCGA
    t=10
    mem=10
    isize=166
    now=$(date +%Y%m%d)

    # targets_bedfile=/mnt/ilustre/users/shuirong.zhang/run/ctdna/config/probe3.hg19.pade0.new.bed
    targets_bedfile=/mnt/ilustre/users/jiasen.wang/danni.li/STR/bed/EGFR.bed
    # factera检测融合基因需要的配置文件
    GRCh37_2bit=/mnt/ilustre/users/yixuezhuanhua/fusion_gene_ann_database/factera/GRCh37.2bit
    GRCh37_bed=/mnt/ilustre/users/yixuezhuanhua/fusion_gene_ann_database/factera/factera-v1.4.4/GRCh37.bed
    # 脚本只支持hg19
    ref=/mnt/ilustre/users/fengbo.zeng/db/genome/human/GRCh37/bwa/GRCh37.fa
    #z-score检测训练组配置文件
    zscore=/mnt/ilustre/users/shuirong.zhang/run/ctdna/config/sample25_mean_std.xls
    #检测样本中snp突变频率(共291个检测样本)
    all_valid_snp_tab_freq=/mnt/ilustre/users/shuirong.zhang/run/ctdna/config/all.valid_snp_tab.freq.xls
}

config(){
    fq1=$s"_R1.fastq.gz"
    fq2=$s"_R2.fastq.gz"
    if [ -f $fq1 ] && [ -f $fq2 ] && [[ $fq1 != $fq2 ]];then
        fqs=($fq1 $fq2)
    else
        fqs=($fq1)
    fi

    echo $now
    echo "sample name: " $s
    echo "fq1: " $fq1
    echo "bedFile: " $targets_bedfile
    echo "ref: " $ref
    echo "cpu: " $t
}

# fq_clean_config(){
#     pre_seq_num6=$(less $fq1|head -4000 |bioawk -c fastx 'NR<=1000{a[substr($seq,7,6)]++}END{for(i in a){print i"\t"a[i]}}'|sort -nk2|tail -1|cut -f2)
#     pre_seq_num9=$(less $fq1|head -4000 |bioawk -c fastx 'NR<=1000{a[substr($seq,7,9)]++}END{for(i in a){print i"\t"a[i]}}'|sort -nk2|tail -1|cut -f2)
#     echo "pre_seq_num6 is : "$pre_seq_num6
#     echo "pre_seq_num9 is : "$pre_seq_num9

#     if [ $pre_seq_num6 -gt 500 ];then
#         if [ $pre_seq_num9 -gt 500 ];then
#             trim_left=15
#         else
#             trim_left=12
#         fi
#     else
#         trim_left=6
#     fi

#     echo "trim left is : "$trim_left
# }

fq_prep(){
    # fq_clean_config
    # for i in 0 1;do
    #     b=$(($i+1))
    #     [ -f ${fqs[$i]} ] && bioawk -c fastx '{print "@"$name" X6:Z:"substr($seq,1,6)"\n"$seq"\n+\n"$qual}' ${fqs[$i]} |\
    #             seqtk trimfq -b $trim_left - > ${s}.trim.${b}.fq
    # done
    
    cutadapt --no-indels -m 15 -a AGTCAGAGANNNNNNTGGAATTCTCGGGTGCCAAGGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT -o ${s}.cut.1.fq.gz -p ${s}.cut.2.fq.gz ${s}_R1.fastq.gz ${s}_R2.fastq.gz
    python /mnt/ilustre/users/jianming.wu/hongyu.chen/add_sixN.py -g $s
    # gzip -c ${s}.clean.1.fq > ${s}_R1.clean.fq.gz
    # gzip -c ${s}.clean.2.fq > ${s}_R2.clean.fq.gz
    # rm ${s}.clean.1.fq ${s}.clean.2.fq ${s}.trim.1.fq ${s}.trim.2.fq
} 

# fq_prep_alternative(){
# 	cutadapt --no-indels -m 15 -a AGTCAGAGANNNNNNTGGAATTCTCGGGTGCCAAGGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG -o ${s}.cut.1.fq.gz ${s}_R1.fastq.gz
# 	cutadapt --no-indels -m 15 -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT -o ${s}.cut.2.fq.gz ${s}_R2.fastq.gz
#     python /mnt/ilustre/users/jianming.wu/hongyu.chen/add_sixN.py -g $s
# }

fq_qc(){
    [ -f $fq1 ] && fastqc -t ${t} $fq1 && fastqc -t ${t} ${s}_R1.clean.fq.gz
    [ -f $fq2 ] && fastqc -t ${t} $fq2 && fastqc -t ${t} ${s}_R2.clean.fq.gz
}

se2bam(){
    bwa mem -C -R '@RG\tID:'$s'\tSM:'$s'\tPL:illumina\tPU:illumina\tLB:illumina' -t $t  $ref $s"_R1.clean.fq.gz" |
        samtools view -@ $t -Sb -q 30 - |
        samtools sort -T $s -@ $t - > $s.bam &&
        samtools index $s.bam
}

pe2bam(){
    bwa mem -C -R '@RG\tID:'$s'\tSM:'$s'\tPL:illumina\tPU:illumina\tLB:illumina' -t $t  $ref $s"_R1.clean.fq.gz" $s"_R2.clean.fq.gz"|
        samtools view -@ $t -Sb  - | samtools sort -T $s -@ $t - > $s.bam
        # samtools view -h $s.raw.bam | samblaster -a -e -d $s.disc.sam -s $s.split.sam -u $s.umc.fasta -o /dev/null &
        # samtools sort -T $s -@ $t  $s.raw.bam > $s.bam
        samtools index $s.bam
}

fq2bam(){
    if [ -f $fq2 ] && [[ $fq1 != $fq2 ]];then
        pe2bam
    else
        se2bam
    fi
}

bam_hit(){
    bedtools intersect -abam $s.bam -b $targets_bedfile  > ${s}.hit_true.bam
    samtools index ${s}.hit_true.bam
    ln -s ${s}.bam ${s}.hit.bam
    ln -s ${s}.bam.bai ${s}.hit.bam.bai
}

split_hit_bam(){
    echo "The hit number of each chromosome: " $(samtools view $s.hit.bam | awk -F"\t" '{a[$3]++}END{for(i in a)print i"\t"a[i]}' |sort -nk1)
    python /mnt/ilustre/users/shuirong.zhang/work/scripts/gada/gada/bin/split_bam.py -i ${s}.hit.bam -1 ${s}.1.bam -2 ${s}.2.bam -3 ${s}.3.bam -4 ${s}.4.bam -5 ${s}.5.bam -6 ${s}.6.bam -7 ${s}.7.bam -8 ${s}.8.bam -9 ${s}.9.bam -0 ${s}.10.bam
    for i in {1..10};do
        if [ -f ${s}.${i}.bam ];then
            echo ${s}.${i}.bam " is existing !"
            samtools index ${s}.${i}.bam &
        else
            echo ${s}.${i}.bam " is not existing !"
        fi 
    done
    wait
}

bam_multi_process_realign(){
#多进程运行,split_hit_bam的输出作为输入
    python /mnt/ilustre/users/shuirong.zhang/work/scripts/gada/gada/bin/multi_process.py -s ${s}
}

merge_realign_bam(){
    samtools view -H -o ${s}.header.sam ${s}.hit.bam
    all_split_bam=$(ls ${s}.realign.*.bam)
    samtools cat -h ${s}.header.sam -o ${s}.hit.realign.merge.bam ${all_split_bam}
    samtools sort -T ${s} -@ ${t} -o $s.hit.realign.bam ${s}.hit.realign.merge.bam
}

clean_multi_process_realign(){
    for i in {1..10};do
        rm ${s}.${i}.bam ${s}.${i}.bam.bai ${s}.realign.${i}.bam ${s}.realign.${i}.bai ${s}.${i}.intervals
    done
    rm ${s}.hit.realign.merge.bam
}

bam_realign(){
    split_hit_bam
    bam_multi_process_realign
    merge_realign_bam
    clean_multi_process_realign
}

bam_dedup(){
    samtools index $s.hit.realign.bam
    samtools index $s.hit.realign.umi_corrected.spolished.bam
    java -Xmx10G -Djava.io.tmpdir=$tmp -jar $picard MarkDuplicates BARCODE_TAG=X6 VALIDATION_STRINGENCY=LENIENT INPUT=$s.hit.realign.umi_corrected.spolished.bam OUTPUT=${s}.markdup.bam METRICS_FILE=${s}.bam.md.metrics
    samtools view -F 1024 -b -o $s.hit.realign.dedup.bam $s.markdup.bam
    samtools sort -o $s.hit.realign.dedup.sort.bam $s.hit.realign.dedup.bam
    rm $s.hit.realign.dedup.bam
    mv $s.hit.realign.dedup.sort.bam $s.hit.realign.dedup.bam
}

bam_filter(){
    samtools view -@ $t -Sb -q 30 $s.hit.realign.dedup.bam > $s.hit.realign.dedup.filter.bam
}
 
bam_qc(){
    bedtools intersect -abam ${s}.hit.realign.dedup.bam -b $targets_bedfile > ${s}.hit_true.realign.dedup.bam
    echo "n_raw: "$(echo $(less ${s}_R1.fastq.gz |wc -l) / 2|bc) > $s.qc
    echo "n_clean: "$(echo $(less ${s}_R1.clean.fq.gz |wc -l) / 2|bc) >> $s.qc
    echo "n_map: "$(samtools view -q 30 -@ $t -c $s.bam) >> $s.qc
    echo "n_map_dedup: "$(samtools view -q 30 -@ $t -c $s.hit.realign.dedup.bam) >> $s.qc  #######去重矫正后map
    echo "n_hit: "$(samtools view -q 30 -@ $t -c ${s}.hit_true.bam) >> $s.qc
    echo "n_hit_dedup: "$(samtools view -q 30 -@ $t -c ${s}.hit_true.realign.dedup.bam) >> $s.qc
    bam_isize
}

bam_prep(){
    bam_hit
    bam_realign
    python /mnt/ilustre/users/jianming.wu/hongyu.chen/cluster.py -g $s
    bam_dedup
    bam_filter
    ln -s $s.hit.realign.dedup.filter.bam $s.valid.bam
    samtools index $s.valid.bam
}

bam2indel(){
    echo -e "$s.valid.bam\t$isize\t$s" > $s.pindel.config
    pindel -f $ref -i $s.pindel.config -c ALL -T $t --include $targets_bedfile -o $s.pindel
    rm $s.pindel_INV         # 有时检测到的INV长度很长(34M多),导致运行pindel2vcf非常慢或者卡住了
    pindel2vcf --max_internal_repeats 5 -r $ref -R GRCh37 -d $now -P $s.pindel -v $s.pindel.indel.vcf

    samtools mpileup -d 100000 -uvf $ref -r 9:133257520-133257521 $s.valid.bam | bcftools call --multiallelic-caller --keep-alts | grep ^#  > $s.header.vcf
    bcftools concat $s.header.vcf $s.pindel.indel.vcf |
        vt sort - | vt uniq - | vt normalize -r $ref - > $s.indel.vcf
}

bam2snp(){
    for bam in $s.valid.bam $s.hit.realign.bam;do
        samtools mpileup -d 100000 -Q 25 -uvf $ref -l $targets_bedfile  -t DP,AD,ADF,ADR,SP $bam |
                bcftools call --multiallelic-caller --keep-alts --targets-file $targets_bedfile --skip-variants indels |
            vt sort - | vt uniq - | vt normalize -r $ref - > $bam.snp.vcf
    done

    ln -s $s.valid.bam.snp.vcf $s.snp.vcf
    ln -s $s.hit.realign.bam.snp.vcf $s.snp0.vcf
}

bam_isize(){
    samtools view -@ $t -F 0x400 $s.valid.bam |
        bioawk -c sam 'BEGIN{print "Id\tTlen\tNumber"}{a[$tlen]++}END{for(i in a){print '$s'"\t"i"\t"a[i]}}' > $s.isize.xls
}

bam2vcf(){
    bam2snp
    bam2indel
}

vcf_ann(){
    for v in snp indel;do
        table_annovar.pl $s.$v.vcf $annovardb -buildver hg19 -protocol refGene,snp138NonFlagged,dbnsfp30a,exac03nontcga,avsnp147 -operation g,f,f,f,f -vcfinput --nastring .
        variant_effect_predictor.pl -i $s.$v.vcf.hg19_multianno.vcf --cache -o $s.$v.vcf.hg19_multianno.vcf.vep.vcf --assembly GRCh37 --offline --force_overwrite --pubmed --vcf
        ln -s $s.$v.vcf.hg19_multianno.vcf.vep.vcf $s.$v.ann.vcf
    done
}

vcf_filter(){
    for v in snp indel;do
        less $s.$v.ann.vcf |
            bio-vcf --sfilter '[s.ad].flatten.min > 3 and r.info["Func.refGene"]=="exonic"' > $s.$v.ann.filter.vcf
        ln -s $s.$v.ann.filter.vcf $s.$v.valid.vcf
    done
}

vcf_prep(){
    vcf_ann
    vcf_filter
}

vcf2tab(){
    for f in $s.snp.vcf $s.snp0.vcf;do
        cat $f |
            bio-vcf --skip-header --eval '[r.chrom,r.pos,r.ref,r.alt.join(","),r.info.dp,r.info.dp4[0..1].reduce(:+),r.info.dp4[2..3].reduce(:+)]' |
            awk -v s=$s '{print s"\t"$0}'> $f.tab

        cat $f |
            bio-vcf --seval 'k=[r.ref, r.alt].flatten.zip([s.ad].flatten).to_h; [k.fetch("A",0), k.fetch("G",0), k.fetch("C",0), k.fetch("T",0)]' |
            awk -v s=$s '{print s"\t"$0}'> $f.tab2
    done

#    for j in $(cat /mnt/ilustre/users/jianming.wu/bin/bed/lvxiaodong.EKB.pos.txt);do grep $j ${s}.snp.vcf.tab ; done > ${s}.result.tab
#	for k in $(cat /mnt/ilustre/users/jianming.wu/bin/bed/lvxiaodong.EKB.pos.txt);do grep $k ${s}.snp0.vcf.tab ; done > ${s}.result.snp0.tab
	cp -f ${s}.snp.vcf.tab ${s}.result.tab
    cp -f ${s}.snp0.vcf.tab ${s}.result.snp0.tab

	awk  -vOFS="\t" '$9=$8+$7 {print}'  ${s}.snp.vcf.tab>${s}.dep.snp.vcf.tab
	cp ${s}.result*.tab /mnt/ilustre/users/jianming.wu/result/ct_str/

    cat $s.indel.ann.vcf |
        bio-vcf --skip-header --set-header "chrom,pos,ref,alt,svtype,func,gene,aa,homseq,ref.dp,alt.dp"  \
                -e  '[chrom,pos,ref,[alt].flatten.join(","), r.info["SVTYPE"], [r.info["Func.refGene"]].flatten.join(","), [r.info["Gene.refGene"]].flatten.join(", "), [r.info["AAChange.refGene"]].flatten.join(", "), r.info.homseq||"."]' \
                --seval 's.ad' |
        awk -v s=$s 'NR==1{print "id\t"$0;next}{print s"\t"$0}'> $s.indel.ann.vcf.tab
}

vcf_qc(){
    echo todo
    echo vcf qc
}

tab_qc(){
    echo "snp_average_depth: "$(awk '{a=a+$7+$8}END{print a/NR}' ${s}.result.snp0.tab) >>$s.qc
    echo "snp_average_depth_dedup: "$(awk '{a=a+$7+$8}END{print a/NR}' ${s}.result.tab) >>$s.qc
    #将质控数据加到数据库dev中
# python /mnt/ilustre/users/shuirong.zhang/work/scripts/gada/gada/bin/add_result2db.py ${s}
}

qc(){
    fq_qc
    bam_qc
    vcf_qc
    tab_qc
}

process_cnv(){
    mkdir -p ${s}.CNV; cd ${s}.CNV; ln -s ../${fq1} ${s}.fastq.gz
    /mnt/ilustre/users/fengbo.zeng/.rvm/rubies/ruby-2.3.0/bin/ruby /mnt/ilustre/users/fengbo.zeng/majorbio/bin/gada.rb cnv -b ~fengbo.zeng/majorbio/data/nchr.20k.gmn.noctdna_panel.bed ${s}
    sed -i '1,6d' ${s}.cnv.sh; sh ./${s}.cnv.sh
    cd ..
}

fusion_gene(){
    mkdir -p ${s}_fusion; cd ${s}_fusion; ln -s ../$s.bam ./; ln -s ../$s.bam.bai ./
    source ~shuirong.zhang/.bashrc.bak
    factera.pl -v -C $s.bam $GRCh37_bed $GRCh37_2bit
    cp -f $s.factera.fusions.txt /mnt/ilustre/users/shuirong.zhang/run/ctdna/2016/result/$s.factera.fusions.xls
    source ~shuirong.zhang/.bashrc
    cd ..
}

valid_snp_anno(){
    /mnt/ilustre/users/shuirong.zhang/work/scripts/gada/gada/bin/ctdna_snp_anno.sh ${s}
    #indel增加一列
    mv ${s}.indel.ann.vcf.tab ${s}.indel.ann.vcf.old.tab
    awk -F"\t" -v OFS="\t" 'NR==1{print $0,"alt.freq"}NR!=1{print $0,$12/($11+$12)}' ${s}.indel.ann.vcf.old.tab > ${s}.indel.ann.vcf.tab
    cp ./${s}.indel.ann.vcf.tab /mnt/ilustre/users/shuirong.zhang/run/ctdna/2016/result/
    rm ${s}.indel.ann.vcf.old.tab
}

table_plot(){
    # 先判断样本是组织还是血液
    sample_type=`grep "sample_name" ${s}.qc |sed 's/sample_name.*FA.*\([A-Z]\{2\}\)/\1/g'`
    if [[ $sample_type == "XZ" ]] || [[ $sample_type == "SL" ]] || [[ $sample_type == "XS" ]]; then    # 组织样本
        s_type="tissue"
    else
        s_type="cfDNA"
    fi
    
    mkdir ${s}_result

    #计算每个基因的覆盖度
    python /mnt/ilustre/users/shuirong.zhang/work/scripts/gada/gada/bin/cal_each_gene_cov.py -t ${s}.result.tab -o ${s}
    cp -f -t ${s}_result ${s}_gene_cov.xls
    #从数据库中提取信息
    python /mnt/ilustre/users/shuirong.zhang/work/scripts/gada/gada/bin/database.py -i ${s}.valid.snp.tab -o ${s}
    #cp -f -t ${s}_result ${s}_C1_pos_ann_result.xls ${s}_C2b_pos_ann_result.xls ${s}_C3_pos_ann_result.xls ${s}_mcgenome_pos_ann.xls ${s}_PGKB_pos_ann.xls

    python /mnt/ilustre/users/shuirong.zhang/work/scripts/gada/gada/bin/ctDNA_report.py $s ctDNA ${s_type}
}

process_result(){
    # 增加snp的突变频率
    awk -F"\t" 'NR==FNR{a[$1"\t"$2"\t"$3"\t"$4]=$5"\t"$6"\t"$7}NR>FNR && FNR==1{print "num\ttotal\tfreq\t"$0}NR>FNR && FNR!=1{if(a[$1"\t"$2"\t"$3"\t"$4]!=""){print a[$1"\t"$2"\t"$3"\t"$4]"\t"$0}else{print "null\tnull\tnull\t"$0}}' ${all_valid_snp_tab_freq} /mnt/ilustre/users/shuirong.zhang/run/ctdna/2016/result/${s}.valid.snp.tab > /mnt/ilustre/users/shuirong.zhang/run/ctdna/2016/result/${s}.valid.snp.tmp.tab
    rm /mnt/ilustre/users/shuirong.zhang/run/ctdna/2016/result/${s}.valid.snp.tab
    mv /mnt/ilustre/users/shuirong.zhang/run/ctdna/2016/result/${s}.valid.snp.tmp.tab /mnt/ilustre/users/shuirong.zhang/run/ctdna/2016/result/${s}.valid.snp.tab

    # 函数目的:将结果重新命名
    guest_name=$(sed -n '8p' ${s}.qc |cut --delimiter=" " -f2)
    sample_name=$(sed -n '9p' ${s}.qc |cut --delimiter=" " -f2)
    cp -R -f ${s}_result /mnt/ilustre/users/shuirong.zhang/run/ctdna/2016/result/
    project_type=`grep "project_type" ${s}.qc|sed 's/project_type.*FA\([0-9]\+\).*/\1/g'`
    # 判断客户名称是否存在qc文件中
    if [ -z $guest_name ];then
        echo "guest name is not existing !"
    else
        cd /mnt/ilustre/users/shuirong.zhang/run/ctdna/2016/result/

        mv ${s}_result ${s}-${sample_name}-${guest_name}_research
        mv $s.factera.fusions.xls ${s}-${sample_name}-${guest_name}.factera.fusions.xls
        mv ${s}.valid.snp.tab ${s}-${sample_name}-${guest_name}.valid.snp.tab
        mv ${s}.indel.ann.vcf.tab ${s}-${sample_name}-${guest_name}.indel.ann.vcf.tab
        mv ${s}.z_score.xls ${s}-${sample_name}-${guest_name}.z_score.xls

        if [ $project_type == 76 ];then
            gene=76
        elif [ $project_type == 8 ];then
            gene=8
        else
            gene=133
        fi
        # 自动对snp和indel进行筛选
        python /mnt/ilustre/users/shuirong.zhang/work/scripts/gada/gada/ctdna/variant_filter.py -g $gene -s ${s}-${sample_name}-${guest_name}.valid.snp.tab -i ${s}-${sample_name}-${guest_name}.indel.ann.vcf.tab
        #文件压缩打包
        zip -r ${s}-${sample_name}-${guest_name}.result.zip ${s}-${sample_name}-${guest_name}_research ${s}-${sample_name}-${guest_name}.factera.fusions.xls ${s}-${sample_name}-${guest_name}.valid.snp.tab ${s}-${sample_name}-${guest_name}.indel.ann.vcf.tab ${s}-${sample_name}-${guest_name}.z_score.xls ${s}-${sample_name}-${guest_name}.valid.snp.filter${gene}.tab ${s}-${sample_name}-${guest_name}.indel.ann.vcf.filter${gene}.tab
        
        if [ $project_type == 8 ];then
            cp -t /mnt/ilustre/users/zhen.yue/result/8gene/ ${s}-${sample_name}-${guest_name}.result.zip
        else
            cp -t /mnt/ilustre/users/zhen.yue/result/other/ ${s}-${sample_name}-${guest_name}.result.zip
        fi

        cd -
        ln -f -s /mnt/ilustre/users/shuirong.zhang/run/ctdna/2016/result/${s}-${sample_name}-${guest_name}.valid.snp.tab ${s}.valid.snp.tab
    fi
    # 研发ctDNA录入到YFA_YCZL表中
    python /mnt/ilustre/users/jiasen.wang/danni.li/STR/programe/addResultTodb_YFA_YCZL_all_valid.py ${s} $targets_bedfile
    python /mnt/ilustre/users/jianming.wu/tmp/time/ctdna-uploaddate.py -i $s

}

NEW_dedup(){
    samtools view $s.bam  > $s.sam
#    /mnt/ilustre/app/pub/perl/5.10.1/bin/perl /mnt/ilustre/users/kefei.huang/program/CTDNA_dedup/dedup_version6/dedup_version6.pl -sam $s.sam -pid_num 10
#    /mnt/ilustre/app/pub/perl/5.10.1/bin/perl /mnt/ilustre/users/kefei.huang/program/CTDNA_dedup/dedup_version6/convert_snp.pl -dedup_vcf $s.call_snp.bam.vcf -raw_vcf $s.snp.vcf                # 修改: $s.snp0.vcf ---> $s.snp.vcf
#   cp annot_result_final_add_ref.xls $s.annot_result_final_add_ref.xls
#    cp annot_result_final.xls $s.annot_result_final.xls
    /mnt/ilustre/app/pub/perl/5.10.1/bin/perl /mnt/ilustre/users/kefei.huang/program/CTDNA_dedup/dedup_version6/dedup_version6.pl -sam $s.sam -sample_name $s -pid_num 10
    /mnt/ilustre/app/pub/perl/5.10.1/bin/perl /mnt/ilustre/users/kefei.huang/program/CTDNA_dedup/dedup_version6/convert_snp.pl -dedup_vcf ${s}_call_snp.bam.vcf -raw_vcf $s.snp.valid.vcf -sample_name $s
    /mnt/ilustre/app/pub/perl/5.10.1/bin/perl /mnt/ilustre/users/kefei.huang/program/CTDNA_dedup/dedup_version6/pindel_new.pl -bam $s.bam -sample_name $s
    /mnt/ilustre/app/pub/perl/5.10.1/bin/perl /mnt/ilustre/users/kefei.huang/program/CTDNA_dedup/dedup_version6/pindel_filter.pl -vcf ${s}_pindel.vcf
}

go(){
    fq_prep
    # fq_prep_alternative
    fq2bam

    bam_prep
    bam2vcf

    vcf_prep
    vcf2tab
    valid_snp_anno

    qc
    fusion_gene
    process_cnv
    table_plot
    process_result
    NEW_dedup
}

clean(){
    rm *fq
    rm *fq.gz
}

if [ -z "$1" ] || [[ $1 == -h ]] || [[ $1 == --help ]]; then
    usage
    exit
elif [[ $1 == "clean" ]];then
    clean
else
    source ~shuirong.zhang/.bashrc
    eval $*
    default
    tmp=$s.tmp && mkdir -p $tmp
    config
    go
fi

